/*******************************************************************************
*
*   Copyright ©2009 Alexandru Marin
*
*   This file is part of nwave/pMRI C library.
*
*   This project started as an MSc degree project during a
*   stage at the Paris-Est university, under joint supervision of prof.
*   Bogdan Ionescu (Image Analysis and Processing Laboratory, "Politehica"
*   University of Bucharest), Caroline Chaux and Jean-Christophe Pesquet
*   (Laboratoire d'Informatique, Université Paris-Est, Institut Gaspard Monge).
*
*   nwave/pMRI C libray is free software: you can redistribute it and/or modify
*   it under the terms of the GNU General Public License as published by
*   the Free Software Foundation, either version 2 of the License, or
*   (at your option) any later version.
*
*   nwave/pMRI C libray is distributed in the hope that it will be useful,
*   but WITHOUT ANY WARRANTY; without even the implied warranty of
*   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*   GNU General Public License for more details.
*
*   You should have received a copy of the GNU General Public License
*   along with nwave/pMRI C library.  If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/

#include "dwt_utils.h"
#include "mri_utils.h"

char* DWT_FILTER_PATH = NULL;

void dwt_filter_set_path(char* path) {

	if (DWT_FILTER_PATH != NULL)
		free(DWT_FILTER_PATH);

	DWT_FILTER_PATH = (char*) malloc((strlen(path) + 1) * sizeof(char));
	assert(DWT_FILTER_PATH!=NULL);

	strcpy(DWT_FILTER_PATH, path);
	printf("* MRI: DWT filter search path set to: %s\n", DWT_FILTER_PATH);
}

dwt_filter_t* dwt_filter_load(char* filtername) {

	if (DWT_FILTER_PATH == NULL)
		dwt_filter_set_path(".");

	char* fullfiltername = mri_strcat(3, DWT_FILTER_PATH, filtername, ".csv");

	FILE* f = fopen(fullfiltername, "r");
	assert(f!=NULL);

	dwt_filter_t* result = dwt_filter_read(f);
	result->name = (char*) malloc((strlen(filtername) + 1) * sizeof(char));
	assert(result->name!=NULL);
	strcpy(result->name, filtername);

	free(fullfiltername);
	fclose(f);

	// Try to read cdjv filters too (in case nothing is found dwt_cdjv_load returns NULL)
	result->cdjv = dwt_cdjv_load(filtername);

	return result;

}

dwt_filter_t* dwt_filter_read(FILE* f) {

	int i;

	dwt_filter_t* result = (dwt_filter_t*) malloc(sizeof(dwt_filter_t));
	assert(result!=NULL);

	// Temp pointer for reading all the filter coefficients at once
	float* temp;
	mri_matrix_descriptor_t desc;
	temp = mri_read_real_matrix(f, &desc);
	assert(desc.ndims == 2);

	result->length = desc.dims[0];

	// Number of rows must be even, because the matrix contains both decomposition and
	// reconstruction filters
	assert(desc.dims[1]%2 == 0);
	result->nbands = desc.dims[1] / 2;

	result->d_filters = (float**) malloc(result->nbands * sizeof(float*));
	assert(result->d_filters!=NULL);
	for (i = 0; i < result->nbands; i++)
		result->d_filters[i] = temp + i * result->length;

	result->r_filters = (float**) malloc(result->nbands * sizeof(float*));
	assert(result->r_filters!=NULL);
	for (i = 0; i < result->nbands; i++)
		result->r_filters[i] = temp + (i + result->nbands) * result->length;

	return result;
}

dwt_filter_t* dwt_filter_alloc(int nbands, int length, char* name) {

	int i;

	dwt_filter_t* result = (dwt_filter_t*) malloc(sizeof(dwt_filter_t));
	assert(result!=NULL);

	result->nbands = nbands;
	result->length = length;

	if (name != NULL) {
		result->name = (char*) malloc(strlen(name) * sizeof(char));
		assert(result->name!=NULL);
		strcpy(result->name, name);
	} else {
		result->name = NULL;
	}

	// Temp pointer for allocating all memory at once
	float* temp = (float*) malloc(2 * nbands * length * sizeof(float));
	assert(temp!=NULL);

	result->d_filters = (float**) malloc(result->nbands * sizeof(float*));
	assert(result->d_filters!=NULL);
	for (i = 0; i < result->nbands; i++)
		result->d_filters[i] = temp + i * result->length;

	result->r_filters = (float**) malloc(result->nbands * sizeof(float*));
	assert(result->r_filters!=NULL);
	for (i = 0; i < result->nbands; i++)
		result->r_filters[i] = temp + (i + result->nbands) * result->length;


	return result;
}

void dwt_filter_write(FILE* f, dwt_filter_t* filter) {

	int i, j;

	fprintf(f, "Filter name: %s\n", filter->name);
	fprintf(f, "Filter Length: %d\n", filter->length);
	fprintf(f, "Number of bands: %d\n", filter->nbands);

	for (i = 0; i < filter->nbands; i++) {

		fprintf(f, "Band #%d decomposition coefficients: [", i + 1);

		// Note that filter coefficients are stored in reverse order
		for (j = filter->length-1; j > 0; j--)
			fprintf(f, "%.2f,", filter->d_filters[i][j]);
		fprintf(f, "%.2f]\n", filter->d_filters[i][j]);
	}

	for (i = 0; i < filter->nbands; i++) {

		fprintf(f, "Band #%d reconstruction coefficients: [", i + 1);

		// Note that filter coefficients are stored in reverse order
		for (j = filter->length-1; j > 0; j--)
			fprintf(f, "%.2f,", filter->r_filters[i][j]);
		fprintf(f, "%.2f]\n", filter->r_filters[i][j]);
	}

}

void dwt_filter_free(dwt_filter_t* filter) {
	// Frees the whole memory chunk allocated for this filter
	free(filter->d_filters[0]);
	free(filter->d_filters);
	free(filter->r_filters);
	free(filter->name);
	if(filter->cdjv != NULL)
		dwt_cdjv_free(filter->cdjv);
	free(filter);
}

dwt_params_t* dwt_params_alloc(int ndims) {

	dwt_params_t* result = (dwt_params_t*) malloc(sizeof(dwt_params_t));
	assert(result!=NULL);

	result->ndims = ndims;

	result->border_strategies = (dwt_border_t*) malloc(ndims * sizeof(dwt_border_t));
	assert(result->border_strategies!=NULL);

	result->filters = (dwt_filter_t**) malloc(ndims * sizeof(dwt_filter_t*));
	assert(result->filters!=NULL);

	result->res_levels = (int*) malloc(ndims * sizeof(int));
	assert(result->res_levels!=NULL);

	return result;
}

void dwt_params_write(FILE* f, dwt_params_t* params) {

	int i;

	fprintf(f, "Number of dimensions: %d\n\n", params->ndims);

	for (i = 0; i < params->ndims; i++) {
		fprintf(f, "Dimension #%d\n", i + 1);
		fprintf(f, "Resolution levels: %d\n", params->res_levels[i]);
		fprintf(f, "Border strategy: %s\n", DWT_BORDER_TO_STR(params->border_strategies[i]));
		dwt_filter_write(f, params->filters[i]);
		fprintf(f, "\n");
	}
}

void dwt_params_free(dwt_params_t* params) {
	free(params->filters);
	free(params->border_strategies);
	free(params->res_levels);
}

void dwt_cdjv_free(dwt_cdjv_t* cdjv) {

	if (cdjv->leftprecond != NULL) {
		free(cdjv->leftprecond[0]);
		free(cdjv->leftprecond);
	}

	if (cdjv->leftprecond != NULL) {
		free(cdjv->rightprecond[0]);
		free(cdjv->rightprecond);
	}

	if (cdjv->leftprecond != NULL) {
		free(cdjv->leftpostcond[0]);
		free(cdjv->leftpostcond);
	}

	if (cdjv->leftprecond != NULL) {
		free(cdjv->rightpostcond[0]);
		free(cdjv->rightpostcond);
	}

	if (cdjv->leftprecond != NULL) {
		free(cdjv->lefthipass[0]);
		free(cdjv->lefthipass);
	}

	if (cdjv->leftprecond != NULL) {
		free(cdjv->leftlopass[0]);
		free(cdjv->leftlopass);
	}

	if (cdjv->leftprecond != NULL) {
		free(cdjv->righthipass[0]);
		free(cdjv->righthipass);
	}

	if (cdjv->leftprecond != NULL) {
		free(cdjv->rightlopass[0]);
		free(cdjv->rightlopass);
	}
}

dwt_cdjv_t* dwt_cdjv_load(char* filtername) {


	FILE* f;
	char* tempname;
	mri_matrix_descriptor_t desc;
	int i, degree;

	// Set default filter path in case not set
	if (DWT_FILTER_PATH == NULL)
		dwt_filter_set_path(".");

	// Read left edge high pass filter - this determines the degree too
	tempname = mri_strcat(3, DWT_FILTER_PATH, filtername, "_lefthipass.csv");
	f = fopen(tempname, "r");
	free(tempname);
	if(f==NULL) {
		return NULL;
	}

	// OK, CDJV exist; allocate result
	dwt_cdjv_t* result = (dwt_cdjv_t*)malloc(sizeof(dwt_cdjv_t));
	assert(result!=NULL);

	float* temp = mri_read_real_matrix(f, &desc);
	fclose(f);

	result->degree = degree = desc.dims[1];
	result->lefthipass = (float**)malloc(degree*sizeof(float*));
	result->leftlopass = (float**)malloc(degree*sizeof(float*));
	result->righthipass = (float**)malloc(degree*sizeof(float*));
	result->rightlopass = (float**)malloc(degree*sizeof(float*));

	result->lefthipass[0] = temp;
	for(i=1;i<degree;i++)
		result->lefthipass[i] = result->lefthipass[0] + i*(3*degree-1);


	// Read left edge lopass filter
	tempname = mri_strcat(3, DWT_FILTER_PATH, filtername, "_leftlopass.csv");
	f = fopen(tempname, "r");
	free(tempname);
	assert(f!=NULL);

	result->leftlopass[0] = mri_read_real_matrix(f, &desc);
	fclose(f);

	for(i=1;i<degree;i++)
		result->leftlopass[i] = result->leftlopass[0] + i*(3*degree-1);


	// Read left edge hipass filter
	tempname = mri_strcat(3, DWT_FILTER_PATH, filtername, "_righthipass.csv");
	f = fopen(tempname, "r");
	free(tempname);
	assert(f!=NULL);

	result->righthipass[0] = mri_read_real_matrix(f, &desc);
	fclose(f);

	for(i=1;i<degree;i++)
		result->righthipass[i] = result->righthipass[0] + i*(3*degree-1);


	// Read left edge lopass filter
	tempname = mri_strcat(3, DWT_FILTER_PATH, filtername, "_rightlopass.csv");
	f = fopen(tempname, "r");
	free(tempname);
	assert(f!=NULL);

	result->rightlopass[0] = mri_read_real_matrix(f, &desc);
	fclose(f);

	for(i=1;i<degree;i++)
		result->rightlopass[i] = result->rightlopass[0] + i*(3*degree-1);




	// Read preconditioning filter
	tempname = mri_strcat(3, DWT_FILTER_PATH, filtername, "_leftprecond.csv");
	f = fopen(tempname, "r");
	free(tempname);

	if(f==NULL) {
		// pre/post cond is not strictly required; in case it's not there
		// just return what has been collected by now
		result->leftprecond = NULL;
		result->rightprecond = NULL;
		result->leftpostcond = NULL;
		result->rightpostcond = NULL;
		return result;
	}

	// OK, we have pre/post conditions
	result->leftprecond = (float**)malloc(degree*sizeof(float*));
	result->rightprecond = (float**)malloc(degree*sizeof(float*));
	result->leftpostcond = (float**)malloc(degree*sizeof(float*));
	result->rightpostcond = (float**)malloc(degree*sizeof(float*));

	result->leftprecond[0] = mri_read_real_matrix(f, &desc);
	fclose(f);

	for(i=1;i<degree;i++)
		result->leftprecond[i] = result->leftprecond[0] + i*degree;


	tempname = mri_strcat(3, DWT_FILTER_PATH, filtername, "_rightprecond.csv");
	f = fopen(tempname, "r");
	free(tempname);
	assert(f!=NULL);

	result->rightprecond[0] = mri_read_real_matrix(f, &desc);
	fclose(f);

	for(i=1;i<degree;i++)
		result->rightprecond[i] = result->rightprecond[0] + i*degree;


	// Read postconditioning filter
	tempname = mri_strcat(3, DWT_FILTER_PATH, filtername, "_leftpostcond.csv");
	f = fopen(tempname, "r");
	free(tempname);
	assert(f!=NULL);

	result->leftpostcond[0] = mri_read_real_matrix(f, &desc);
	fclose(f);

	for(i=1;i<degree;i++)
		result->leftpostcond[i] = result->leftpostcond[0] + i*degree;


	tempname = mri_strcat(3, DWT_FILTER_PATH, filtername, "_rightpostcond.csv");
	f = fopen(tempname, "r");
	free(tempname);
	assert(f!=NULL);

	result->rightpostcond[0] = mri_read_real_matrix(f, &desc);
	fclose(f);

	for(i=1;i<degree;i++)
		result->rightpostcond[i] = result->rightpostcond[0] + i*degree;



	return result;
}

void dwt_cdjv_write(FILE* f, dwt_cdjv_t* cdjv) {

	int i, j;

	fprintf(f, "\nCohen, Daubechies, Jawerth and Vial border wavelets\n");

	fprintf(f, "Degree: %d\n", cdjv->degree);

	fprintf(f, "\nPrecondition left matrix:\n");
	for(i=0;i<cdjv->degree;i++) {
		for(j=0;j<cdjv->degree;j++)
			fprintf(f, " %.3f ", cdjv->leftprecond[i][j]);
		fprintf(f, ";");
	}
	fprintf(f, "\nPrecondition right matrix:\n");
	for(i=0;i<cdjv->degree;i++) {
		for(j=0;j<cdjv->degree;j++)
			fprintf(f, " %.3f ", cdjv->rightprecond[i][j]);
		fprintf(f, ";");
	}

	fprintf(f, "\nLeft hipass matrix:\n");
	for(i=0;i<cdjv->degree;i++) {
		for(j=0;j<3*cdjv->degree-1;j++)
			fprintf(f, " %.3f ", cdjv->lefthipass[i][j]);
		fprintf(f, ";");
	}
	fprintf(f, "\nLeft lopass matrix:\n");
	for(i=0;i<cdjv->degree;i++) {
		for(j=0;j<3*cdjv->degree-1;j++)
			fprintf(f, " %.3f ", cdjv->leftlopass[i][j]);
		fprintf(f, ";");
	}

	fprintf(f, "\nRight hipass matrix:\n");
	for(i=0;i<cdjv->degree;i++) {
		for(j=0;j<3*cdjv->degree-1;j++)
			fprintf(f, " %.3f ", cdjv->righthipass[i][j]);
		fprintf(f, ";");
	}

	fprintf(f, "\nRight lopass matrix:\n");
	for(i=0;i<cdjv->degree;i++) {
		for(j=0;j<3*cdjv->degree-1;j++)
			fprintf(f, " %.3f ", cdjv->rightlopass[i][j]);
		fprintf(f, ";");
	}

	fprintf(f, "\nPostcondition left matrix:\n");
	for(i=0;i<cdjv->degree;i++) {
		for(j=0;j<cdjv->degree;j++)
			fprintf(f, " %.3f ", cdjv->leftpostcond[i][j]);
		fprintf(f, ";");
	}
	fprintf(f, "\nPostcondition right matrix:\n");
	for(i=0;i<cdjv->degree;i++) {
		for(j=0;j<cdjv->degree;j++)
			fprintf(f, " %.3f ", cdjv->rightpostcond[i][j]);
		fprintf(f, ";");
	}


}
